1  Aires protégées

1.1 Données de l’association Vahatra

Les études sur les aires protégées s’appuient fréquemment sur la base WDPA (World Database on Protected Area), consultable en ligne sur https://protectedplanet.net. On s’aperçoit dans le cas de Madagascar que cette base de données comporte de nombreuses erreurs (qu’on étudiera plus bas). La base rassemblée par l’association Vahatra dans le cadre de la monographie qu’elle a coordonnée sur l’ensemble des aires protégées terrestres malgaches semble beaucoup plus fiable (Goodman et al. 2018). Les données en question sont disponibles sur le portail https://protectedareas.mg avec une licence creative commons (CC-BY).

Le bloc de code ci-dessous (cliquer sur “code” pour visualiser), présente la séquence d’opérations réalisées pour préparer les données.Pour comprendre certaines opérations contenues dans le bloc de code, il est utile d’être familier de la syntaxe de R et des packages du tidyverse. Voir le chapitre Chapter 10.

Code

library(tidyverse)
library(lubridate)
library(sf)
library(tmap)
library(geodata)
library(cowplot)
library(wdpar)
library(gt) # Pour faciliter le rendu des tableaux (et ils sont jolis)

# Le shapefile est composé d'une série de fichiers, (.shp, .dbf, .prj, .shx)
# qui doivent avoir le même nom et être au même endroits pour être ouverts en
# même temps. Comme souvent, ils sont compressés ensemble dans un fichier zip.
# On commence par dézipper (décompresser) ce fichier.
unzip("data/Vahatra98AP.zip", exdir = "data/Vahatra")
# On importe dans R en pointant vers le fichier .shp, mais c'est bien toute la
# collection de fichiers homonymes .shp, .dbf, .shx qui est chargée.
AP_Vahatra <- st_read("data/Vahatra/Vahatra98AP.shp", quiet = TRUE) %>%
  # Il manque la projection (pas de fichier .prj), on la spécifie à la main
  st_set_crs("EPSG:4326") # EPSG 4325 = WSG 84 = le standard pour le web

# L'option ci-dessous est un peu cryptique : des caractéristiques topologiques
# de la carte source sont incompatibles avec la possibilité d'avoir des objets
# sphériques dans sf. Cela disparait si on désactive cette possibilité
sf_use_s2(FALSE) 

# Identification des dates ----------------------------------------------------
# Cette section est un brin complexe, à base de manipulation de chaînes de 
# caractères et de dates

# Détecte les dates écrites 2 avril 2020 ou 02 avril 2020, etc.
date_ecrite <- "[:digit:]{1,2} [:alpha:]* [:digit:]{4}"
# Détecte les dates écrites 02/04/20 ou 02.04.20 ou 02.04.2020, etc.
date_abrev <- "[:digit:]{2}[:punct:][:digit:]{2}[:punct:][:digit:]{2,4}"
# Des années écrites à 2 chiffres
date_ecrite_an_abrev <- "[:digit:]{1,2} [:alpha:]* [:digit:]{4}"
# Détecte l'une ou l'autre des formes précédentes
toute_date <- paste(date_ecrite, date_abrev, date_ecrite_an_abrev, sep = "|")
# Détecte une mention d'année seule : 1984, 2015, etc.
annee_seule <- "[:digit:]{4}"
# Détecte les formes indicatrices d'un changement
mention_changement <- "Changement|changement|anciennement|actuel|auparavant"
# Une fonction qui traduit les dates écrites en toutes lettre du français à 
# l'anglais (pour les parser ensuite car ça ne fonctionne qu'en anglais)
trad_dates <- function(date_fr) {
  str_replace_all(date_fr,
                  c("janvier" = "January",
                    "fevrier" = "February",
                    "mars" = "March",
                    "avril" = "April",
                    "mai" = "May",
                    "juin" = "June",
                    "juillet" = "July",
                    "aout" = "August",
                    "septembre|setembre" = "September",
                    "octobre" = "October",
                    "novembre" = "November",
                    "decembre|decmbre" = "December"))
}
# Cette fonction remplace 01.04.58 par 01.04.1958 et marche avec . ou /
# On indique avec limite le nombre d'année où on considère que c'est 1900 vs 2000
complete_annee <- function(date_abrev, limite = 20) {
  if (str_detect(date_abrev, "([:punct:])([:digit:]{2})[:punct:]?$")) {
    date_abrev <- str_remove(date_abrev, ":punct:]?$")
    if (as.numeric(str_extract(date_abrev, "[:digit:]{2}$")) > limite) {
      date_abrev <- str_replace(date_abrev, 
                                "([:punct:])([:digit:]{2})[:punct:]?$", "\\119\\2")
    } else {
      date_abrev <- str_replace(date_abrev, 
                                "([:punct:])([:digit:]{2})[:punct:]?$", "\\120\\2")
    }
  }
  return(date_abrev)
}
# La fonction précédente est unitaire, on la transforme pour qu'elle s'applique à une liste.
complete_liste_dates <- function(liste_dates) {
  map(liste_dates, complete_annee)
}

AP_Vahatra <- AP_Vahatra %>%
  # On extrait les dates des champs de texte
  mutate(date_creation = str_extract_all(creation, toute_date), 
         # Une date a un format incohérent, on la recode à la main
         date_creation = ifelse(creation == "Créée le 07 aout 04",
                                "07 aout 2004", date_creation),
         date_creationA = map(date_creation, 1), # La 1ère date
         date_creationB = map(date_creation, 2)) %>% # Si 2 dates, la seconde
  # On traduit les mois en anglais pour une conversion au format date
  mutate(across(c("date_creationA", "date_creationB"), trad_dates)) %>%
  mutate(across(c("date_creationA", "date_creationB"), complete_liste_dates)) %>%
  mutate(across(c("date_creationA", "date_creationB"), dmy)) %>%
  mutate(date_creation = case_when(is.na(date_creationB) ~ date_creationA,
                                    date_creationA > date_creationB ~ date_creationB,
                                    date_creationA <= date_creationB ~ date_creationA),
         date_modification = case_when(is.na(date_creationB) ~ date_creationB,
                                       date_creationA < date_creationB ~ date_creationB,
                                       date_creationA >= date_creationB ~ date_creationA),
         # On repère si il y a eu un changement de statut ou de frontières
         mention_changement = str_detect(creation, mention_changement)) %>%
    # On enlève les colonnes inutiles
  select(-date_creationA, -date_creationB) %>%
  # On place les colonnes créées à gauche pour les inspecter facilement
  relocate(date_creation:mention_changement, .after = creation) 

# Après une vérification manuelle, on remarque les données de l'association Vahatra comportent des mentions incomplètes pour certaines aires, qui n'ont pas été extraites:


# Lokobe : 31 décembre 1927
AP_Vahatra <- AP_Vahatra %>%
  mutate(date_creation = case_when(nom == "Lokobe" ~ ymd("1927-12-31"),
                                   nom == "Mantadia" ~ ymd("1989-01-11"),
                                   TRUE ~ date_creation),
         date_modification = case_when(nom == "Bemaraha" ~ ymd("2011-07-06"),
                                       nom == "Lokobe" ~ ymd("2011-07-06"),
                                       nom == "Tsaratanana" ~ ymd("2011-07-06"),
                                       nom == "Pic d'Ivohibe" ~ ymd("2015-04-28"),
                                       nom == "Mantadia" ~ ymd("2002-08-07"),
                                       TRUE ~ date_modification)) %>%
  st_make_valid() # fiabilise qu'il n'y a pas d'erreurs topologiques
# dir.create("AP_Vahatra")
# st_write(AP_Vahatra, "out/AP_Vahatra.shp")
# writexl::write_xlsx(st_drop_geometry(AP_Vahatra), "AP_Vahatra.xlsx")

Le bloc de code suivant génère une carte interactive. On a également inclus des lignes de code qui permettent de formater la carte joliment pour un rendu figé (pdf/LaTeX, html statique, word), mais ce code est “commenté”, c’est-à-dire qu’on a placé des dièses au début de chaque ligne, de sorte qu’il ne s’exécute pas (R n’exécute jamais ce qui se trouve à droite d’un # sur une ligne). Pour plus de détails sur la manière dont on produit des cartes, voire l’annexe : Cartes simples en R

Code

if (file.exists("data/contour_mada.rds")) {
  load("data/contour_mada.rds")
} else {
  contour_mada <- gadm(country = "Madagascar", resolution = 1, level = 0,
                     path = "data/GADM") %>%
  st_as_sf()
# On enregistre contour_mada pour s'en servir par la suite
save(contour_mada, file = "data/contour_mada.rds")
}

# On génère un rendu cartographique
tmap_mode("view") # En mode interactif
tm_shape(contour_mada) +
  tm_borders() +
  tm_shape(AP_Vahatra) + 
  tm_polygons(col = "cat__iucn", alpha = 0.6, title = "Catégorie IUCN") +
  tmap_options(check.and.fix = TRUE) # +

Code

# Les dièses en début de ligne font que ce qui suit ne s'exécute pas.
# La suite est uniquement pour les rendus fixes (tmap_mode = "plot"), p. ex. pour les pdf
  # # NB : on note les positions en majuscules quand on veut coller aux marges
  # tm_credits("Sources: WDPA et GADM", position = c("RIGHT", "BOTTOM"),
  #            size = 0.6) +
  # tm_layout(main.title = "Aires protégées de Madagascar",
  #           # NB : position en minuscules pour laisser un espace avec la marge
  #           main.title.position = c("center", "top"),
  #           main.title.size = taille_titres_cartes,
  #           legend.position = c("left", "top"),
  #           legend.outside = TRUE)

On peut également réaliser un graphique qui présente l’historique de création des aires protégées. Pour plus de précisions sur la manière de produire des graphiques en R, voir l’annexe correspondante.

Code

# On ordonne les nom d'aires protégées dans l'ordre de leur séquence de création
ordre_chrono_AP <- AP_Vahatra %>%
  arrange(desc(date_creation), desc(nom)) %>%
  pull(nom)
# On transforme le champ "nom" de caractère, à une catégorisation ordonnée où
# l'ordre correspond 
AP_Vahatra_carte <- AP_Vahatra %>%
  mutate(nom = factor(nom, levels = ordre_chrono_AP),
         cat_taille = case_when(hectares > 300000 ~ 2,
                                hectares > 150000 ~ 1.5,
                                hectares >  50000 ~ 1,
                                             TRUE ~ 0.5)) %>%
  rename(`Catégorie IUCN` = cat__iucn)

# On crée un graph pour les anciennetés
graph_gauche <- AP_Vahatra_carte %>%
  ggplot(aes(x = date_creation, xend = ymd("2022-10-01"), y = nom, yend = nom, 
                   color = `Catégorie IUCN`)) +
  geom_segment(size = 2) +
  ggtitle("Ancienneté") +
  theme(axis.ticks.y = element_blank(),
        legend.position = "none",
        axis.title = element_blank(),
        axis.text.x.bottom =  element_text(angle = 45, hjust = 1),
        axis.text.x.top = element_text(angle = -45, hjust = 1)) + 
  scale_x_date(sec.axis = dup_axis())

graph_droite <- AP_Vahatra_carte %>%
  ggplot(aes(x = 0, xend = hectares/100, y = nom, yend = nom, 
                   color = `Catégorie IUCN`)) +
  geom_segment(size = 2) + 
  ggtitle("Surface (km2)") +
  theme(axis.text.y = element_blank(),
        axis.title = element_blank(),
        axis.text.x.bottom =  element_text(angle = 45, hjust = 1),
        axis.text.x.top = element_text(angle = -45, hjust = 0),
        legend.position = "none") + 
  scale_x_continuous(sec.axis = dup_axis())

legende <- get_legend(graph_gauche  +
                        guides(color = guide_legend(nrow = 1)) +
                        theme(legend.position = "bottom"))

# On colle les deux
graphs <- plot_grid(graph_gauche, graph_droite, rel_widths = c(2.2, 1),
          nrow = 1)
plot_grid(graphs, legende, ncol = 1,
          rel_heights = c(1,.1))

Il faut aussi s’assurer qu’on filtre bien les entités analysées selon un critère pertinent. Actuellement, on exclut les aires marines. Il pourrait toutefois sembler utile d’écarter les aires dont le statut de protection est considéré comme trop faible. Il pourrait aussi être pertinent de ne garder que les aires protégées comportant un niveau minimum de couvert forestier : autrement, cela signifie que la forêt n’est pas un habitat pertinent pour les écosystèmes que la démarche de conservation cherche à protéger dans cette aire.

1.2 World Database on Protected Areas

On commence par télécharger et présenter ces données.

Code

# On regarde si les données WDPA sont disponibles sur l'ordinateur qui exécute
if (file.exists("data/WDPA/WDPA_Oct2022_MDG-shapefile.zip")) {
  # Si oui, on charge
  WDPA_Mada <- wdpa_read("data/WDPA/WDPA_Oct2022_MDG-shapefile.zip")
} else {
  # Si non, on télécharge depuis protectedplanet
  WDPA_Mada <- wdpa_fetch("Madagascar", wait = TRUE,
                      download_dir = "data/WDPA") 
}



# TODO: inclure graph et description des données.

1.3 Comparaison des données Vahatra et WDPA

On commence par visualiser les différences spatiales entre les polygones, en affichant les 10 qui sont les plus différents entre les WDPA et Vahatra.

Code

# On harmonise les noms qui sont parfois notés différemment entre les sources
AP_Vahatra <- AP_Vahatra %>%
  mutate(nom_wdpa = case_when(
    nom == "Corridor Forestier Bongolava" ~ "Corridor forestier Bongolava",
    nom == "Ranobe PK32" ~ "Ranobe PK 32",
    str_detect(nom, "Ambositra-Vondrozo") ~ "Corridor Forestier Ambositra Vondrozo",
    nom == "Réserve deTampolo" ~ "Réserve de Tampolo",
    nom == "Bombetoka Beloboka" ~ "Bombetoka Belemboka",
    nom == "Ampananganandehibe-Behasina" ~ "Ampanganandehibe-Behasina",
    nom == "Forêt Sacrée Alandraza Analavelo" ~ "Analavelona", # vérfié sur carte : les mêmes
    nom == "Réserve speciale Pointe à Larrée" ~ "Réserve spéciale Pointe à Larrée", 
    nom == "Vohidava-Betsimalaho" ~ "Vohidava Betsimalao", 
    nom == "Anjanaharibe Sud" ~ "Anjanaharibe_sud",
    nom == "Iles Radama/Sahamalaza" ~ "Sahamalaza Iles Radama",
    nom == "Kalambatritra" ~ "Kalambatrika",
    nom == "Mananara-Nord" ~ "Mananara Nord",
    nom == "Kirindy - Mitea" ~ "Kirindy Mite",
    nom == "Midongy du Sud" ~ "Befotaka Midongy", # Vérifié sur la carte
    nom == "Montagne d'Ambre/Forêt d'Ambre" ~ "Montagne d'Ambre",
    nom == "Tsimanampesotsa" ~ "Tsimanampesotse",
    nom == "Pic d'Ivohibe" ~ "Ivohibe",
    nom == "Forêt Naturelle de Petriky" ~ "Forêt Naturel de Petriky",
    nom == "Tsingy de Namoroka" ~ "Namoroka",
    nom == "Réserve de Ressources Naturelle Mahimborondro" ~ "Mahimborondro",
    str_detect(nom, "Complexe Tsimembo Manambolomaty") ~ "Complexe Tsimembo Manambolomaty",
    nom == "Mandrozo" ~ "Zone Humide de Mandrozo",
    nom == "Paysage Harmonieux Protégés Bemanevika" ~ "Complexe des Zones Humides de Bemanevika",
    nom == "Nord Ifotaky" ~ "INord fotaky",
    TRUE ~ nom)) %>%
  arrange(nom_wdpa) %>%
  mutate(rownum = row_number())

# On sauvegarde le résultat
save(AP_Vahatra, file = "data/ch1_AP_Vahatra.rds")

# On ne garde que les aires de WDPA qui apparaissent dans Vahatra
WDPA_commun <- WDPA_Mada %>%
  filter(NAME %in% AP_Vahatra$nom_wdpa) %>%
  filter(!(NAME == "Analalava" & IUCN_CAT == "Not Reported")) %>%
  filter(!(NAME == "Site Bioculturel d'Antrema" & IUCN_CAT == "Not Reported")) %>%
  filter(DESIG != "UNESCO-MAB Biosphere Reserve") %>%
  arrange(NAME)  %>%
  mutate(rownum = row_number())
       
# Cette fonction calcule la part d'un polygone incluse dans un 
# autre polygone et retourne un ratio entre 0 et 1
ratio_inclus <- function(x, y) {
  inclus <- st_intersection(x, y)
  ratio <- st_area(inclus) / st_area(x)
  return(ratio)
}

# On calcule la part des polygones Vahatra incluse dans les polgones WDPA 
V_in_W <- map2_dbl(WDPA_commun$geometry, AP_Vahatra$geometry, ratio_inclus)
# Puis l'inverse
W_in_V <- map2_dbl(AP_Vahatra$geometry, WDPA_commun$geometry, ratio_inclus)
# On fait un facteur des deux
recoupement_mutuel <- V_in_W * W_in_V
# Qu'on ramène dans les jeux de données d'origine
WDPA_commun2 <- bind_cols(WDPA_commun, V_in_W = V_in_W, W_in_V = W_in_V,
                         recoupement_mutuel = recoupement_mutuel) %>%
  arrange(recoupement_mutuel, rownum)
AP_Vahatra2 <- bind_cols(AP_Vahatra, V_in_W = V_in_W, W_in_V = W_in_V,
                        recoupement_mutuel = recoupement_mutuel) %>%
  arrange(recoupement_mutuel, rownum)

# On prend maintenant les 5 les plus éloignés et on les visualise
min_recoup <- WDPA_commun2 %>%
  filter(row_number() <= 10) %>%
  select(nom_wdpa = NAME, rownum) %>%
  mutate(source = "WDPA") %>%
  bind_rows(select(filter(AP_Vahatra2, rownum %in% .$rownum), nom_wdpa, rownum)) %>%
  mutate(source = ifelse(is.na(source), "Vahatra", source))
tmap_mode("plot")
min_recoup %>%
  tm_shape() +
  tm_polygons() +
  tm_facets(by = c("nom_wdpa", "source")) +
  tm_layout(panel.label.size=3)

On peut également comparer ceux pour lesquels on a des différences de date ou de statut.

Code

# On garde seulement les métadonnées qu'on veut comparer
WDPA_a_comparer <- WDPA_commun %>% # On repart des AP communes
  st_drop_geometry() %>% # Plus besoin de spatial
  select(nom_wdpa = NAME, type_wdpa = INT_CRIT, cat_iucn_wdpa = IUCN_CAT,
         year_wdpa = STATUS_YR) # On ne garde que les colonnes à comparer

verif_meta_wdpa <-AP_Vahatra %>%
  st_drop_geometry() %>% # Pas besoin d'un jeu spatial
  select(nom:date_modification, nom_wdpa) %>% # colonnes à garder dans Vahatra
  # On renomme la catégorie IUCN de Vahatra et on code les NA comme dans WDPA
  mutate(cat_iucn = ifelse(is.na(cat__iucn), "Not Reported", cat__iucn)) %>%
  relocate(cat_iucn, .before = cat__iucn) %>% # Nouvelle colonne près de l'ancienne
  left_join(WDPA_a_comparer, by = "nom_wdpa") %>% # On rassemble Vahatra et WDPA
  select(-nom_wdpa, -cat__iucn) %>% # On enlève les colonnes inutiles
  # On compare les dates et statuts
  mutate(`Différence de date` = year(date_creation) != year_wdpa,
         `Différence de statut` = cat_iucn != cat_iucn_wdpa)

verif_meta_wdpa %>%
  summarise(`Nombre d'aires protégées comparées` = n(),
            `Différence de date` = sum(`Différence de date`),
            `Différence de statut` = sum(`Différence de statut`)) %>%
  gt() %>%
  tab_header(title = paste("Différences entre les données de WDPA et celles de",
                     "l'assciation Vahatra sur les aires protégées terrestres",
                     "à Madagascar"))

Différences entre les données de WDPA et celles de l'assciation Vahatra sur les aires protégées terrestres à Madagascar
Nombre d'aires protégées comparées	Différence de date	Différence de statut
98	68	40

Dans les cas qu’on peut comparer, les données de l’association Vahatra semblent plus fiables. On va donc privilégier l’utilisation de ces dernières.

On va également visualiser les aires de WDPA qui ne sont pas contenues dans Vahatra.

Code

WDPA_exclu <- WDPA_Mada %>%
  filter(!(NAME %in% AP_Vahatra$nom_wdpa))

tmap_mode("view")
WDPA_exclu %>%
  tm_shape() +
  tm_polygons(col = "IUCN_CAT")

Code

ratio_terrestre <- function(x) {
  inclus <- st_intersection(x, contour_mada$geometry)
  ratio <- st_area(inclus) / st_area(x)
  return(ratio)
}

# On crée un grand polygone avec toutes les AP dans Vahatra
AP_Vahatra_fusion <- st_union(AP_Vahatra)

# On calcule pour les aires protégées de WDPA qui ne sont pas dans Vahatra
# dans quelle mesure elles sont terrestres et pas superposées à d'autres AP
# déjà dans Vahatra
WDPA_exclu <- WDPA_exclu %>%
  filter(st_geometry_type(.) == "MULTIPOLYGON") %>%
  mutate(part_terrestre = map2_dbl(.$geometry, contour_mada$geometry, 
                                   ratio_inclus),
         part_deja_autre = map2_dbl(.$geometry, AP_Vahatra_fusion, 
                                    ratio_inclus))

# On garde celles qui sont au moins 25% terrestre et 75% pas superposées
WDPA_a_inclure <- WDPA_exclu %>%
  filter(part_terrestre >= 0.25 & part_deja_autre <= 0.25) %>%
  mutate(full_name = paste(INT_CRIT, NAME)) %>%
  select(nom = NAME, WDPAID, full_name, creation = STATUS_YR)

On a visiblement des aires protégées qu’il serait pertinent d’inclure et qui ne sont pas dans Vahatra.

1.4 Enjeux de fiabilité des données d’aires protégées

Important pour l’analyse : si périmètres pas juste => phénomènes de leakage, faux positifs ou faux négatifs.

Enjeu aussi des métadonnées : date ou type sont importants pour l’analyse et celle-ci perd en fiabilité si ces informations ne sont pas correctes.

2 Caractéristiques spatiales

Elle finit termine par enregistrer le jeu de données produit sur la machine qui exécute le code. Ce bloc ne s’exécute que si le jeu de données résultant n’est pas détecté sur la machine. Si le jeu de données résultant du script précédent est déjà disponible sur la machine, alors le bloc précédent ne s’exécute pas et celui qui s’exécute est le suivant

library(tidyverse)
library(mapme.biodiversity)
library(sf)

if (file.exists("data/Vahatra_poly.rds")) {
  load("data/Vahatra_poly.rds")
} else {

  Vahatra_poly <- AP_Vahatra %>%
    filter(st_geometry_type(.) == "MULTIPOLYGON") %>%
    st_cast("POLYGON")
  
  # Constitution d'un portefeuille (voir la documentation)
  Vahatra_poly <- init_portfolio(x = Vahatra_poly, 
                                 years = 2000:2020,
                                 outdir = "data/mapme_Vahatra",
                                 cores = 24,
                                 add_resources = TRUE,
                                 verbose = TRUE)
  
  # Données d'accessibilité de Nelson et al. (2018)
  Vahatra_poly <-  get_resources(x = Vahatra_poly, resource = "nelson_et_al",  
                                 range_traveltime = "5k_110mio")
  # Modèle numérique de terrain SRTM de la NASA
  Vahatra_poly <- get_resources(x = Vahatra_poly , resource = "nasa_srtm")
  
    # Indicateurs d'accessibilité
  Vahatra_poly <- calc_indicators(x = Vahatra_poly,
                                  "traveltime",  stats_accessibility = "mean",
                                  engine = "extract")
  # Indicateurs de relief de terrain
  Vahatra_poly <- calc_indicators(x = Vahatra_poly,
                                  indicators = c("tri", "elevation"),
                                  stats_tri = "mean", stats_elevation = "mean")
  
  #   # On récupère aussi les données de Global Forest Watch sur le couver forestier
  # Vahatra_poly <- get_resources(x = Vahatra_poly, 
  #                               resources = c("gfw_treecover", "gfw_lossyear"))
  #   # Indicateurs de couvert forestier
  # Vahatra_poly <- calc_indicators(x = Vahatra_poly,
  #                                 indicators = "treecover_area", 
  #                                 min_cover = 30, min_size = 1)
  
  save(Vahatra_poly, file = "data/Vahatra_poly.rds")
}

Mapme produit des colonnes imbriquées pour chaque observation, car dans bien des cas, on peut avoir plusieurs valeurs (par année) pour une même observation, voire plusieurs variables (par exemple, le calcul de l’indicateur traveltime produit des estimations de distance par rapport à une ville pour plusieurs tailles de ville possible. Lorsqu’on spécifie une taille, il produit deux variables : la distance estimée et la taille de la ville prise en compte pour l’estimation.

Cette imbrication n’est pas indispensable pour les trois variables calculées ici (indice de terrain accidenté, distance à une ville et altitude), car on ne cherche qu’une valeur par observation. On va donc dés-imbriquer les variables.

# Valeur agrégées par AP (moyennes pondérées par la surface)
Vahatra_vars_terrain <- Vahatra_poly %>%
  unnest(cols = c(tri, elevation, traveltime)) %>%
  st_drop_geometry() %>%
  select(nom, hectares, indice_accidente = tri_mean, dist_ville = minutes_mean, 
         altitude = elevation_mean) %>%
  group_by(nom) %>%
  summarise(indice_accidente = weighted.mean(indice_accidente, hectares,
                                             na.rm = TRUE),
            dist_ville = weighted.mean(dist_ville, hectares,
                                       na.rm = TRUE),
            altitude = weighted.mean(altitude, hectares,
                                     na.rm = TRUE))
# Valeurs qu'on insère dans le jeu de données de travail
load("data/ch1_AP_Vahatra.rds")
AP_Vahatra <- AP_Vahatra %>%
  left_join(Vahatra_vars_terrain, by = "nom")

save(AP_Vahatra, file = "data/ch2_AP_Vahatra.rds")

On doit aussi se rappeler que les aires protégées sont parfois composées de plusieurs polygones disjoints et que mapme.biodiversity a calculé chaque indicateur pour chaque polygone séparément. Pour chaque aire protégée, on va donc faire la moyenne de ces indicateurs, pondérée par la surface respective de chaque polygone.

Données d’accessibilité : attention car elles présentent un possible biais d’endogénéité. La construction de route au cours des dernières décennies peut être lié à l’établissement ou non d’aires protégées. L’inclusion d’une variable de contrôle qui peut être en partie affectée par notre variable de traitement (la conservation) est susceptible de problème. Il existe une carte de 2000 qui pourrait être mobilisée :

On notera que plusieurs autres indicateurs peuvent être calculés à partir du pabkage mapme.biodiversity:

• active_fire_counts: Calculate active fire counts based on NASA FIRMS polygonsactive_fire_properties: Calculate active fire properties based on NASA FIRMS polygons

• biome: Calculate biomes statistics (TEOW) based on WWF

• drought_indicator: Calculate drought indicator statistics

• ecoregion: Calculate terrestrial ecoregions statistics (TEOW) based on WWF

• landcover: Calculate area of different landcover classes

• mangroves_area: Calculate mangrove extent based on Global Mangrove Watch (GMW)

• population_count: Calculate population count statistics (Worldpop)

• precipitation_chirps: Calculate precipitation statistics based on CHIRPS

• precipitation_wc: Calculate precipitation statistics

• soilproperties: Calculate Zonal Soil Properties

• temperature_max_wc: Calculate maximum temperature statistics

• temperature_min_wc: Calculate minimum temperature statistics based on WorldClim

• traveltime: Calculate accessibility statistics

• treecover_area: Calculate treecover statistics

• treecover_area_and_emissions: Calculate treeloss statistics

• treecoverloss_emissions: Calculate emission statistics

• tri: Calculate Terrain Ruggedness Index (TRI) statistics

3 Couvert forestier

3.1 Carvalho et al. (source MODIS)

Pour commencer, on récupère le travail réalisé par Carvalho et al. (2020) qui complète les informations physiques de Goodman et al. (2018) avec des données relatives au couvert forestier en 1996, 2006 et 2016 et la diversité d’espèces présentes.

library(tidyverse)
library(readxl)
library(sf)

# Voir le chapitre "Fondamentaux R" pour une aide à l'import.
sup2 <- read_xlsx("data/Carvalho2018sup2.xlsx", # Enplacement du fichier
                  skip = 8, # Premières lignes du tableau excel à ne pas lire
                  n_max = 101,  # on ne lit pas les dernières lignes (notes)
                  col_types = c("text", "text", "text", "text", 
        "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", "numeric"))
sup4 <- read_xlsx("data/Carvalho2018sup4.xlsx", skip = 6,
                  col_types = c("text", "numeric", "text", "numeric", "numeric", 
        "numeric", "numeric", "numeric", "numeric", "numeric"))
        

# Carvalho et al. 2008 document in their supp. material 2: "The three parcels that made up
# Andohahela (Parcels I, II and III) comprised different types of dominant vegetation and
# associated animal species, and were exposed to distinct pressures. Andohahela was analysed
# in its entirety (site number 57), as well as separated"

sup2 <- sup2 %>% 
  mutate(PA = recode(PA, `Andohahela complete` = "Andohahela"),
         num_atlas_ = as.integer(`Site number`))

sup4 <- sup4 %>%
  filter(`Habitat type` == "TOTAL") %>%
  mutate(num_atlas_ = as.numeric(Parcel))


load("data/ch2_AP_Vahatra.rds")

AP_Vahatra <- AP_Vahatra %>%
  left_join(sup2, by = "num_atlas_") %>%
  relocate(PA, .after = nom) %>%
  left_join(sup4, by = "num_atlas_")

On complète cette information avec des données de couvert forestier.

3.2 Mapme (exemple GFC)

La procédure de traitement de ces fichiers sur Mapme est analogue à celle employée dans la section Chapter 2.

#| eval: false

# On charge les polygones travaillés au chapitre 2
load("data/Vahatra_poly.rds")

# Constitution d'un portefeuille (voir la documentation)
AP_poly <- init_portfolio(x = AP_poly, 
                          years = 2000:2020,
                          outdir = "data_s3/mapme",
                          cores = 16,
                          add_resources = TRUE,
                          verbose = TRUE)


# Get GFW data
AP_poly <- get_resources(x = AP_poly, 
                         resources = c("gfw_treecover", "gfw_lossyear", 
                                       "gfw_emissions"))

# Données d'accessibilité de Nelson et al. (2018)
AP_poly <-  get_resources(x = AP_poly, resource = "nelson_et_al",  
                          range_traveltime = "5k_110mio")
# Données de qualité des sols (uniquement teneur )
AP_poly  <-  get_resources(x = AP_poly,
                           resources = "soilgrids",  layers = "clay", 
                           depths = "5-15cm", stats = "mean")
# Modèle numérique de terrain SRTM de la NASA
AP_poly <- get_resources(x = AP_poly, resource = "nasa_srtm")
# Données de feux
AP_poly <- get_resources(x = AP_poly, resource = "nasa_firms",
                         instrument = "MODIS")

# Get 
# Indicateurs de couvert forestier
AP_poly  <- calc_indicators(x = AP_poly,
                            indicators = "treecover_area_and_emissions", 
                            min_cover = 30, min_size = 1)

# On réorganise les variables -----------------------------------------------
# Couvert forestier annuel en colonnes plutôt qu'en liste imbriquée  
deforest_par_an <- AP_poly %>%
  unnest(treecover_area_and_emissions) %>%
  select(assetid, nom, years, treecover) %>%
  filter(!is.na(years)) %>%
  pivot_wider(names_from = "years", values_from = "treecover", 
              names_prefix = "treecover_") %>%
  st_drop_geometry()

# Jointure du couvert forestier en colonne avec les données initiales
AP_poly <- AP_poly %>%
  select(-treecover_area_and_emissions) %>%
  left_join(deforest_par_an, by = "assetid")

Toutefois, en raison d’un problème liés à la gestion des calculs volumineux, les calculs pour certaines aires protégées renvoient des données aberrantes. Ce point sera mis à jour dans ce guide dès la résolution des erreurs rencontrées.

A ce stade on se concentrera donc sur les données de TFM présentées plus bas.

3.3 Google Earth Engine (exemple GFC)

La plateforme Google Earth Engine est un outil particulièrement pratique et performant pour mobiliser et traiter des données satellitaires. Google Earth Engine peut être utilisé :

• en interrogeant son API, et notamment :
  • en python, avec la librairie gee permet d’interroger l’API de Google Earth Engine.
  • en R, au travers de la librairie rgee. Cette dernière est relativement facile d’usage, mais elle est difficile à configurer. Pour aller plus loin : https://r-earthengine.com/rgeebook/
• directement sur la plateforme https://code.earthengine.google.com/

La consolde de codage de Google Earth Engine prend la forme suivante :

Diagramme des composants de la console Google Earth Engine

Le langage utilisé sur cet interface est du Javascript. Ci-dessous un exemple de code qui génère les surface (en hectares) de perte de couvert forestier. Pour fonctionner, ce code doit être collé dans un script sur la plateforme Google Earth Engine lancé en cliquant sur “Run”, puis en cliquant sur “Tasks” pour exécuter le code.

#| eval: false

 scale = 30
    
    // PREPARE DATA
    //look at tree cover, find the area
    var treeCover = gfc2021.select(['treecover2000']);
    var areaCover = treeCover.multiply(ee.Image.pixelArea())
                    .divide(10000).select([0],["areacover"])
    // total loss area
    var loss = gfc2021.select(['loss']);
    var areaLoss = loss.gt(0).multiply(ee.Image.pixelArea())
                   .divide(10000).select([0],["arealoss"]);
    // total gain area
    var gain = gfc2021.select(['gain'])
    var areaGain = gain.gt(0).multiply(ee.Image.pixelArea())
                   .divide(10000).select([0],["areagain"]);
    // final image
    var total = gfc2021.addBands(areaCover)
                .addBands(areaLoss)
                .addBands(areaGain)

    // TOTAL COVER
    // Map cover area per feature
    var districtSums = areaCover.reduceRegions({
      collection: testgu,
      reducer: ee.Reducer.sum(),
      scale: scale,
    });         
    
    var addVar = function(feature) {

      // function to iterate over the sequence of years
      var addVarYear = function(year, feat) {
        // cast var
        year = ee.Number(year).toInt()
        feat = ee.Feature(feat)

        // actual year to write as property
        var actual_year = ee.Number(2000).add(year)

        // filter year:
        // 1st: get mask
        var filtered = total.select("lossyear").eq(year)
        // 2nd: apply mask
        filtered = total.updateMask(filtered)

        // reduce variables over the feature
        var reduc = filtered.reduceRegion({
          geometry: feature.geometry(),
          reducer: ee.Reducer.sum(),
          scale: scale,
          maxPixels: 1e9
        })

        // get results
        var loss = ee.Number(reduc.get("arealoss"))
        var gain = ee.Number(reduc.get("areagain"))

        // set names
        var nameloss = ee.String("loss_").cat(actual_year)
        var namegain = ee.String("gain_").cat(actual_year)

        // set properties to the feature
        return feat.set(nameloss, loss, namegain, gain)
      }

      // iterate over the sequence
      var newfeat = ee.Feature(years.iterate(addVarYear, feature));

      // return feature with new properties
      return newfeat
    }

    // Map over the FeatureCollection
    var areas = districtSums.map(addVar);
    
    // Export PA deforestation to a CSV file.
    Export.table.toDrive({
      collection: areas,
      description: 'forest_loss_WDPA_Madagascar',
      fileFormat: 'CSV'
    });

3.4 Python (exemple TMF)

Fichiers préparés en python (code à venir), directement sur les rasters.

#| eval: false
# On charge les fichiers préparés par Marc en python (contient 2 feuilles)
tableur_tmf <- "data/TMFchangeYear_AP_Vahatra.xlsx"
# On commence par charger la feuille déforestation
tmf_vahatra_defor <- read_excel(tableur_tmf,
                              sheet = "TMFdeforestationYear") %>%
  select(nom, starts_with("TMF")) # on ne garde que les variables pertinentes
# On fait ensuite de même avec la feuille dégradation
tmf_vahatra_degrad <- read_excel(tableur_tmf,
                              sheet = "TMFdegradationYear") %>%
  select(nom, starts_with("TMF")) # Onn ne garde que les feuilles pertinentes

AP_Vahatra <- AP_Vahatra %>%
  left_join(tmf_vahatra_degrad, by = "nom")
  # left_join(tmf_vahatra_defor, by = "nom") %>%


TMF_ratios <- AP_Vahatra %>%
  st_drop_geometry() %>%
  select(nom, starts_with("Forest cover"), starts_with("TMF")) %>%
  pivot_longer(cols = starts_with("TMF"), 
               names_to = "variable", 
               values_to = "surface_ha") %>%
  mutate(an_valeur = str_extract(variable, "[:digit:]{4}"),
         an_valeur = as.numeric(an_valeur),
         surface_ratio = case_when(
           an_valeur < 2000 ~ surface_ha / `Forest cover (ha) in 2006`,
           an_valeur > 2009 ~ surface_ha / `Forest cover (ha) in 2016`,
           TRUE ~ surface_ha / `Forest cover (ha) in 2006`),
         variable = str_replace(variable, "HA", "ratio")) %>%
  select(nom, variable, surface_ratio) %>%
  pivot_wider(names_from = variable, values_from = surface_ratio) %>%
  select(nom, starts_with("TMF"))
  
AP_Vahatra <- AP_Vahatra %>%
   left_join(TMF_ratios, by = "nom")

save(AP_Vahatra, file = "data/ch3_AP_Vahatra.rds")

3.5 Alternatives

Si on n’est pas à l’aise avec les outils mentionnés plus haut, l’outil Geoquery d’AidData permet d’obtenir des statistiques par aire administrative. Il est également possible de formuler des demandes spécifiques pour d’autres polygones que des aires administratives au travers d’un formulaire dédié.

4 Méthode randomisée

ATTENTION : Il va de soi que les AP malgaches n’ont à aucun moment été assignées aléatoirement. Lors de cette séquence, on fait “comme si” pour montrer la manière dont les données sont analysées quand il y a eu assignation aléatoire. On verra en fin de session les limites d’une telle approche et dans les suivantes des manières de construire des contrefactuels plus vraisemblables pour un sujet comme celui-ci.

Le jeu de données AP_Vahatra contenait 98 aires protégées, avec des géométries de types “multi-polygones”. Certaines de ces aires protégées étaient en effet composées de plusieurs polygones disjoints. Ces polygones disjoints ont été scindés pour être traités séparément dans le jeu de données AP_poly. Avant de repasser sur des analyses agrégées, on va agrégéer les statistiques d’AP_poly afin d’avoir pour chaque variable une valeur par aire protégée.

library(tidyverse)
library(gt)
library(broom)
library(stargazer)
library(sf)
library(lubridate)

load("data/ch3_AP_Vahatra.rds")

rct_AP_Mada <- AP_Vahatra %>%
  st_drop_geometry() %>%
  rename(`Déforestation 1996-2006 (%)` = 
           `Forest loss (ha) between 1996-2006 (percent loss)`,
         `Déforestation 2006-2016 (%)` = 
           `Forest loss (ha) between 2006-2016 (percent loss)`,
         `Surface (ha)` = hectares) %>%
  mutate(Groupe = ifelse(year(date_creation) < 2015, "Traitement", "Controle"),
         `Couvert forestier en 1996 (%)` = `Forest cover (ha) in 1996` / 
                                              `Surface (ha)` * 100,
         `Déforestation 1996-2016 (%)` = 
           (`Forest loss (ha) between 1996-2006 (absolute loss)` + 
             `Forest loss (ha) between 2006-2016 (absolute loss)`) /
           `Forest cover (ha) in 1996` * 100)

# On fait une série de tests de comparaison de moyenne
t_tests <- rct_AP_Mada %>% 
  # On applique aux variables de déforestation, couvert en 96 et taille
  summarise(across(ends_with("(%)") | ends_with("(ha)"),# toutes finissent ainsi
                   ~ t.test(.[Groupe == "Controle"], # on applique un t.test
                            .[Groupe == "Traitement"])$p.value)) %>%
  mutate(Groupe = "t-test")

equilibre_avant <- rct_AP_Mada %>%
  group_by(Groupe) %>%
  summarise(`Nombre d'aires` = n(),
            `Sans forêt` = sum(is.na(`Couvert forestier en 1996 (%)`)), 
            `Surface (ha)` = mean(`Surface (ha)`),
            `Couvert forestier en 1996 (%)` = 
              mean(`Couvert forestier en 1996 (%)`, na.rm = TRUE)) %>%
  bind_rows(t_tests) %>% # On colle tous les t-tests 
  mutate(across(!Groupe, ~round(., 2))) %>%# arrondit tout sauf colonne "Groupe"
  select(-starts_with("Déforestation"))# On ne garde que les t-tests de baseline

# Ce qui suit est une série d'opération pour formater le rendu en tableau
equilibre_avant %>%
  t() %>% # On transpose lignes <=> colonnes
  as.data.frame() %>% # La transposition a altéré le format, on remet en tableau
  tibble::rownames_to_column() %>% # On met le nom des lignes en 1° colonne
  # "Truc pour renommer avec le contenu de la première ligne
  `colnames<-` (filter(., row_number() == 1)) %>% 
  filter(row_number() != 1)%>% # Enlève la 1° ligne qui est maintenant en entête
  gt() %>%
  tab_header(title = "Equilibre des variables avant intervention",
             subtitle = "(exercice : \"comme si\" c'était une RCT)") %>%
  tab_source_note("Source : Association Vahatra et Carvalho et al. 2018")

On a à première vue des équilibres limités “avant intervention”. En moyenne, les deux groupes sont assez proches en termes de surface et de couvert forestier et le test de Stydent ne permet pas de rejeter l’hypothèse nulle.

On va maintenant s’intéresser aux différences de déforestation observées “après intervention” dans le groupe de traitement.

comparaison_apres <- rct_AP_Mada %>%
  group_by(Groupe) %>%
  summarise(across(starts_with("Déforestation"), ~ mean(., na.rm = TRUE))) %>%
  bind_rows(t_tests) %>% # On colle tous les t-tests 
  mutate(across(!Groupe, ~round(., 2))) %>%# arrondit tout sauf colonne "Groupe"
  select(Groupe, starts_with("Déforestation"))# On ne garde que les t-tests de baseline


# Même procédure que plus haut pour formater le rendu en tableau
comparaison_apres  %>%
  t() %>% # On transpose lignes <=> colonnes
  as.data.frame() %>% # La transposition a altéré le format, on remet en tableau
  tibble::rownames_to_column() %>% # On met le nom des lignes en 1° colonne
  # "Truc pour renommer avec le contenu de la première ligne
  `colnames<-` (filter(., row_number() == 1)) %>% 
  filter(row_number() != 1)%>% # Enlève la 1° ligne qui est maintenant en entête
  gt() %>%
  tab_header(title = "Moyennes après intervention",
             subtitle = "(exercice : \"comme si\" c'était une RCT)") %>%
  tab_source_note("Source : Association Vahatra et Carvalho et al. 2018")

On peut également réaliser une régression simple.

def_96_06 <- lm(`Déforestation 1996-2006 (%)`  ~ Groupe, data = rct_AP_Mada)
def_06_16 <- lm(`Déforestation 2006-2016 (%)`  ~ Groupe, data = rct_AP_Mada)
def_96_16 <- lm(`Déforestation 1996-2016 (%)`  ~ Groupe, data = rct_AP_Mada)

tidy(def_96_06) %>%
  gt() %>%
  tab_header(title = "Déforestation 1996-2006 (%)",
             subtitle = "(exercice : \"comme si\" c'était une RCT)") %>%
  tab_source_note("Source : Association Vahatra et Carvalho et al. 2018")

tidy(def_06_16) %>%
  gt() %>%
  tab_header(title = "Déforestation 2006-2016 (%)",
             subtitle = "(exercice : \"comme si\" c'était une RCT)") %>%
  tab_source_note("Source : Association Vahatra et Carvalho et al. 2018")

tidy(def_96_16) %>%
  gt() %>%
  tab_header(title = "Déforestation 1996-2016 (%)",
             subtitle = "(exercice : \"comme si\" c'était une RCT)") %>%
  tab_source_note("Source : Association Vahatra et Carvalho et al. 2018")

On tente une autre mise en forme

stargazer(def_96_06, def_06_16, def_96_16, type = "text",
          title = "Impact de la conservation sur la perte de couvert forestier",
          notes = "Données : Association Vahatra et Carvalho et al. 2018")

References

Carvalho, Fabio, Kerry A. Brown, Adam D. Gordon, Gabriel U. Yesuf, Marie Jeanne Raherilalao, Achille P. Raselimanana, Voahangy Soarimalala, and Steven M. Goodman. 2020. “Methods for Prioritizing Protected Areas Using Individual and Aggregate Rankings.” Environmental Conservation 47 (2): 113–22. https://doi.org/10.1017/S0376892920000090.

Goodman, Steven M., Marie Jeanne Raherilalao, Sébastien Wohlhauser, Jean Clarck N. Rabenandrasana, Herivololona M. Rakotondratsimba, Fanja Andriamialisoa, and Malalarisoa Razafimpahanana. 2018. Les Aires Protégées Terrestres de Madagascar: Leur Histoire, Description Et Biote. Association Vahatra.